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Abstract  The  sensitivity  of  the  bifurcation  of  the  North 
Equatorial  Current  in  the  Pacific  to  different  wind  pro¬ 
ducts  is  investigated.  Variations  of  the  bifurcation  lati¬ 
tude  with  season  is  simulated  in  a  purely  wind-driven 
model  and  is  found  to  be  in  agreement  with  recent 
observations.  The  seasonal  cycle  is  nearly  independent 
of  the  wind  climatology,  but  the  annual  average  latitude 
depends  on  the  wind  stress  curl.  It  is  also  shown  that 
in  the  upper  ocean,  the  poleward  shift  in  bifurcation 
latitude  with  depth  is  realistic  in  our  simple  model.  This 
implies  that  given  a  stratification  close  to  the  observed, 
it  is  primarily  the  wind  forcing  that  determines  the 
location  of  the  bifurcation  and  its  seasonal  variation. 

Keywords  North  Equatorial  Current  •  Bifurcation  • 
North  Pacific  •  Numerical  modeling 


1  Introduction 

The  Pacific  North  Equatorial  Current  (NEC")  bifurcates 
into  the  northward-flowing  Kuroshio  and  southward¬ 
flowing  Mindanao  Current  at  the  Philippine  coast.  An 
important  role  of  these  two  western  boundary  currents 
is  closing  the  mass  budgets  of  the  subtropical  and 
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tropical  gyres,  respectively,  and  the  NEC  bifurcation 
marks  the  separation  of  the  two  gyres  at  the  western 
boundary.  It  is  therefore  of  interest  to  understand  the 
processes  that  control  the  location  of  the  latitude. 


IT  Observations 

The  bifurcation  latitude  varies  seasonally  and  with 
depth.  Toole  et  al.  (1990)  noted  a  bifurcation  latitude 
of  14.1°  in  September  1987,  while  it  was  observed  to 
be  at  12.6°  in  April  1988,  suggesting  a  seasonal  change, 
and  Nitani  (1972)  first  reported  a  poleward  shift  in  bi¬ 
furcation  latitude  with  depth.  Based  on  geostrophic  cal¬ 
culations  relative  to  1,500  dbar  from  the  Levitus  (1982) 
annual  mean  climatology,  Ou  et  al.  (1998)  showed  that 
the  NEC  volume  transport  of  41  Sv  bifurcates  into  a 
northward  transport  of  14  Sv  and  southward  transport 
of  27  Sv:  they  also  confirmed  a  northward  shift  in 
bifurcation  latitude  with  depth  from  about  13.5°  N  at 
the  surface  to  18°  N  at  500  dbar  in  good  agreement  with 
the  water  property  distributions  (Ou  et  al.  1999). 

In  a  recent  comprehensive  reanalysis  of  historical  hy¬ 
drography  data,  Qu  and  Lukas  (2003)  determined  that 
the  bifurcation  for  currents  averaged  over  the  upper 
1,000  m  occurs  at  17.2°  N  in  December  and  at  14.8°  N 
in  July.  Its  annual  mean  location  was  at  15.2°  N,  close 
to  the  estimate  from  Sverdrup  theory  (14.6°  N),  varying 
from  13.3°  N  near  the  surface  to  about  20°  N  at  1,000  m. 

Ou  and  Lukas  (2003)  found  the  bifurcation  latitude 
to  be  closest  to  the  equator  in  June-July  and  at  its 
northernmost  position  in  December  (Fig.  1).  Through¬ 
out  the  year,  the  bifurcation  latitude  increases  with 
depth,  in  particular  during  the  boreal  winter.  Note  that 
the  bifurcation  from  observations  is  well-defined  only 
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Fig-1  Meridional  current 
component  along  the  coast  of 
the  Philippines  in  July  (top) 
and  January  (bottom).  Unit  is 
centimeters  per  second.  The 
bifurcation  latitudes  from  the 
layer  model  forced  by 
ECMWF  winds  are  shown  by 
filled  circles  (adapted  from 
Qu  and  Lukas  2(X)3) 


J  ul 


Jan 


at  depths  less  than  about  500  m  during  the  summer  and 
to  about  600  m  during  the  winter. 

The  motivation  for  this  study  is  the  above-mentioned 
analysis  by  Qu  and  Lukas  (2(X)3),  who  determined  the 
annual  eycle  of  the  bifurcation  based  on  hydrographie 
observations.  Their  results  are  further  confirmed  by 
recent  results  from  a  high-resolution  ocean  general 
circulation  model  (OGCM)  by  Kim  et  al.  (2004),  but 
differed  significantly  from  earlier  modeling  results  by 
Qiu  and  Lukas  (1996),  who  used  a  simple  wind-driven 
reduced  gravity  model  and  suggested  that  the  bifurca¬ 
tion  was  determined  by  the  wind  stress  eurl  alone.  The 
discrepancy  between  the  work  of  Qiu  and  Lukas  (1996) 
and  the  later  research  questions  to  what  extent  a  sim¬ 
ple  wind-driven  model  is  able  to  capture  the  essential 
dynamics  needed  to  reproduce  a  realistic  bifurcation  or 
if  additional  vertical  modes  plays  a  role  or  if  buoyancy 
forcing  is  needed.  It  also  raises  the  question  whether 
the  bifurcation  latitude  is  a  robust  feature  that  is  closely 
related  to  the  gyre  circulation,  or  if  it  is  sensitive  to  local 
winds  and  boundary’  flows. 

In  this  work,  we  address  these  questions.  We  use 
wind-driven  layer  models  forced  by  different  wind  cli¬ 


matologies  and  compare  the  annual  cycle  of  the  bifur 
cation  latitude  with  the  observational  results  from  Qu 
and  Lukas  (2003). 

1 .2  Models 

The  model  studies  by  Qiu  and  Lukas  (1996)  using  a 
linear  reduced  gravity  model  were  forced  by  the  Florida 
State  University  (FSIJ)  wind  stress  (1961-1992)  and 
showed  that  the  northern  extreme  occurred  in  October 
and  the  southernmost  position  occurred  in  February. 
That  integration  was  close  to  being  in  phase  with  the 
annual  march  of  the  latitude  of  zero  wind  stress  curl, 
but  with  much  reduced  amplitude.  Applying  a  nonlin¬ 
ear  numerical  model,  Qiu  and  Lukas  (1996)  found  a 
somewhat  larger  amplitude  than  in  the  linear  model 
and  a  shift  in  phase  so  the  southernmost  latitude  was 
reached  in  Mareh-April  and  the  northernmost  latitude 
in  Oetober-November.  Both  their  linear  and  non-linear 
model  results  differ  from  the  observations  of  Qu  and 
Lukas  (2003)  discussed  in  the  previous  section. 

Qiu  and  Lukas  (1996)  also  found  that  in  the  annual 
mean,  the  bifurcation  latitude  was  shifted  about  0.5°  to 
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the  south  of  the  annual  mean  of  the  zero  wind  stress 
curl  latitude.  From  theory  (McCreary  and  Lu  1994), 
it  is  seen  that  in  case  of  a  linear  model  with  an  active 
subsurface  layer,  the  bifurcation  latitude  will  in  general 
not  be  located  at  the  latitude  of  zero  wind  stress  curl 
even  in  case  of  steady  winds.  They  found  that  the 
bifurcation  latitude  for  the  transport  field  lies  at  the 
zero  wind  stress  curl  line  only  for  an  idealized  wind 
stress  field  that  is  independent  of  latitude.  It  is  therefore 
not  clear  from  theory  at  which  latitude  the  bifurcation 
should  occur  for  flow  in  a  realistic  ocean  model  forced 
by  seasonal  winds. 

1.3  Present  research 

In  this  paper,  we  demonstrate  that  the  observed  change 
in  bifurcation  latitude  with  season  as  well  as  the  pole- 
ward  shift  with  depth  in  the  upper  ocean  can  be  re¬ 
produced  in  a  wind-driven  layer  model.  This  implies 
that  annual  buoyancy  forcing  does  not  play  an  essential 
role  for  the  bifurcation  above  the  main  thermocline. 
On  decadal  time  scales  or  longer,  buoyancy  force  is 
essential  in  maintaining  the  stratification.  It  is  also  of 
interest  to  determine  if  the  bifurcation  is  sensitive  to 
different  wind  climatologies,  so  a  number  of  diverse 
wind  products  were  used,  and  we  also  briefly  explore 
if  the  shape  of  the  coast  line,  local  and  remote  forcing, 
or  bottom  topography  plays  an  important  role. 

2  Numerical  model 

The  ocean  model  is  based  on  the  multi-layer  upper- 
ocean  model  (Jensen  1991,  1993),  extended  to  include 
mixed  layer  physics  and  an  option  to  include  bottom 
topography  (Jensen  1998,  2(X)1).  In  the  models  used 
here,  lateral  stresses  are  proportional  to  the  local  de¬ 
formation  rate  of  the  flow  (Smagorinsky  1963)  rather 
than  constants  as  in  previous  versions.  A  Smagorinsky 
constant  of  0.2  was  used.  Details  about  the  model  is 
given  in  Jensen  (2(X)1). 

The  model  domain  covers  the  global  ocean  north  of 
60°  S  to  60°  N  with  a  horizontal  resolution  of  1/4°  and 
closed  poleward  boundaries.  The  configuration  used 
for  most  experiments  is  a  4.5-layer  reduced  gravity 
model:  It  has  four  active  layers  with  the  average  initial 
thickness  of  80,  120,  250,  and  600  m  for  layers  1  to  4, 
respectively.  The  deepest  layer  is  infinitely  deep  and  is 
at  rest.  The  sigma-theta  levels  are  chosen  to  be  23.6, 
25.4,  26.5,  27.2,  and  28.2  and  are  held  constant  in  each 
isopycnal  layer.  In  order  to  maintain  a  realistic  layer 
thickness,  upper  limits  and  lower  limits  are  imposed 
on  each  layer  before  mass  is  entrained  or  detrained 


from  adjacent  layers.  This  is  the  hydromixing  approach 
implemented  in  other  layer  models  (e.g.,  Hurlburt  et  al. 
1996).  The  layer  thickness  limits  without  entrainment 
are  50-120  m  for  layer  1,  80-200  m  for  layer  2,  1 20— 
250  m  for  layer  3,  and  200-800  m  for  layer  4.  Details 
of  the  method  of  constraining  layer  thickness  as  used  in 
this  model  are  given  in  Jensen  (2001). 

The  reduced  gravity  model  is  the  simplest  confi¬ 
guration  that  extends  the  experiments  by  Qiu  and 
Lukas  (1996)  to  include  a  change  of  the  bifurcation 
with  depth.  The  advantage  of  this  model  is  the  relative 
low  cost  of  integration,  in  spite  of  the  high  horizontal 
resolution  needed  for  a  good  resolution  of  the  bifur¬ 
cation  latitude.  The  relative  simple  model  allows  us 
to  perform  many  experiments  otherwise  prohibitive 
due  to  high  computational  costs.  A  five-layer  model 
configuration  used  to  test  potential  impact  of  bot¬ 
tom  topography  is  described  in  a  subsequent  section. 
Differences  in  the  flows  between  the  4.5-  and  five-layer 
model  configurations  turn  out  to  be  relatively  small  in 
the  upper  fours  layers  in  the  area  of  interest,  so  we  can 
use  the  purely  baroclinic  flows  from  the  4.5-layer  model 
as  reference. 


3  Forcing  for  the  control  run  and  model  spin-up 

Climatological  monthly  mean  wind  stress  from  the 
European  Centre  for  Medium  Range  Forecast 
(ECMWF)  reanalysis  (ERA-15)  is  used  as  forcing  for 
a  control  run.  The  time  period  from  1979  to  1988 
corresponding  to  the  first  Atmospheric  Modeling  Inter¬ 
comparison  Project  wras  chosen.  The  wind  stress  is  in¬ 
terpolated  from  the  2.5°  spatial  resolution  to  the  ocean 
model  grid  using  cubic  splines  and  linear  interpolation 
in  time  so  it  is  updated  each  model  time  step. 

A  complete  spin-up  of  the  global  ocean  requires  in¬ 
tegration  over  several  hundred  years  and  is  not  feasible. 
By  using  a  density  stratification  based  on  the  annual 
average  found  in  the  vicinity  of  the  bifurcation  region 
computed  from  the  World  Ocean  Atlas  94  (Levitus 
et  al.  1994;  Levitus  and  Boyer  1994),  it  is  assumed 
that  a  thermal  adjustment  is  not  needed.  Hence,  there 
is  no  buoyancy  forcing,  so  the  model  is  purely  wind- 
driven.  The  adjustment  time  for  the  model  is  set  by 
propagation  of  Rossby  waves  across  the  basin  (e.g.,  Gill 
1982,  pp.  507-512).  The  adjustment  time  for  vertical 
mode  n  is 

=  Lf2/ficl  (1) 

where  L  is  the  width  of  the  basin,  /  the  Coriolis  pa¬ 
rameter  at  a  given  latitude,  fi  its  variation  with  lati¬ 
tude,  and  c„  the  internal  wave  speed  of  vertical  mode 
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n.  For  the  tropical  and  subtropical  Pacific,  the  dom¬ 
inant  vertical  mode  is  the  first  baroclinic  mode.  ITie 
stratification  above  implies  C\  =  3.6  m/s,  but  using  a 
more  conservative  estimate  of  C\  =  3  m/s,  allowing  for 
layer  adjustment  to  slow  the  waves,  and  a  basin  scale 
of  L  =  15,000  km  gives  a  spin-up  time  of  6.2  years 
at  20°  N  as  a  minimum.  Higher  vertical  modes  have 
much  longer  adjustment  times;  for  instance,  the  second 
internal  mode  wave  speed  is  1.62  m/s,  requiring  a  mini¬ 
mum  adjustment  time  of  21  years.  However,  the  higher 
modes  are  also  damped  by  friction  which  limits  the 
adjustment  time  significantly  for  mode  3  and  higher  and 
prevents  influence  over  large  longitude  ranges.  Figure  2 
shows  the  baroclinic  sea  surface  height  p  computed  for 
a  4.5-layer  model  (e.g.,  Jensen  2001,  2003): 

>1  =  ^— (//,-//«)  (2) 
rr  as 


where  px  is  the  density,  //;  the  layer  thickness,  and  //(>, 
the  initial  layer  thickness  at  rest  of  layer  /.  This  surface 
height  includes  contributions  from  all  four  baroclinic 
modes  and  provides  a  measure  of  the  gyre  circulation 
and  its  magnitude.  As  expected,  the  major  circulation 
features  in  the  tropical  are  present  after  6  years.  The 
subtropical  gyres  expand  as  the  model  continue  to  spin- 
up,  although  the  differences  are  relative  small  between 
year  12  and  year  20  in  the  tropics. 

A  concise  discussion  of  basin  adjustment  times  is 
given  in  McCreary  et  al.  (2007).  Using  the  same  damp¬ 
ing  rate  as  in  their  work,  the  mode-2  damping  time  scale 
is  very  close  to  the  mode-2  adjustment  time  for  Rossby 
waves,  i.e.,  20  years. 

We  use  the  circulation  after  20  years  of  spin-up  as  the 
control  case.  While  this  is  too  short  to  completely  re¬ 
move  spin-up  effects  on  the  ocean  gyres,  it  is  sufficient 
to  reproduce  the  basic  flow  features  in  the  tropical  and 
subtropical  gyres.  The  model  is  eddy  permitting  with  a 


Fig.  2  Sea  surface  height  displacement  from  initial  condition  using  a  3-month  average  centered  on  January  6  on  year  6  {(op),  year  12 
(center),  and  year  20  (bottom)  during  the  spin-up.  Unit  is  meters 
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flow  dependent  low  eddy  viscosity,  so  a  90-day  running 
average  was  made  to  the  model  fields  before  analysis. 


4  Seasonal  flow  and  transports 

For  the  ECMWF,  control  run  transports  of  the  upper 
850  m  of  three  sections  were  computed.  One  section  is 
from  the  Philippine  coast  to  140°  E  along  18°  N.  On 
the  annual  average,  a  northward  Kuroshio  transport 
of  13.5  Sv  is  found.  The  second  section  is  across  140° 
E  between  8°  N  and  18°  N.  The  NEC  transport  in 
this  section  has  an  annual  average  of  36.1  Sv.  The 
third  section,  from  the  coast  along  8°  N  to  140°  E, 
provides  an  estimate  of  the  Mindanao  Current  trans¬ 
port  of  14.1  Sv,  a  magnitude  slightly  above  the  model 
Kuroshio  transport.  Tozuka  et  al.  (2002)  also  re¬ 
ported  a  near  equal  division  of  the  NEC  transport  to 
the  Kuroshio  Current  and  Mindanao  Current.  Esti¬ 
mates  from  observational  studies  have  found  that  the 
Mindanao  transport  exceeds  the  Kuroshio  transport; 
for  instance,  Qu  et  al.  (1998)  reported  it  just  over  twice 
as  large,  while  Kashino  et  al.  (2009)  found  it  larger 
by  29%  and  60%  during  winters  of  2006  and  2008, 
respectively.  The  model  only  includes  the  upper  850  m 
of  the  ocean  in  the  bifurcation  region,  and  there  is  a  net 
subduction  through  the  lowest  layer  of  8.5  Sv. 

Figure  3  shows  the  monthly  transport  from  the  three 
sections.  We  find  that  the  covariation  between  the  NEC 
and  the  Mindanao  Current  is  greater  than  between  the 


NEC  and  Kuroshio  on  a  seasonal  time  scale.  Simula¬ 
tions  with  an  ocean  general  circulation  model  shows 
similar  seasonal  correlations  (Kim  et  al.  2(X)4).  The 
currents  in  the  layer  model  simulation  agree  well  with 
the  observed  flow  (Qu  et  al.  1998,  1999)  and  flow  in 
general  circulation  models  (Kim  et  al.  2(X)4).  Figure  4 
shows  the  monthly  mean  currents  when  the  bifurcation 
is  in  its  northernmost  position  in  January  and  in  its 
southernmost  position  in  July.  These  results  increase 
our  confidence  in  the  simple  wind-driven  model  for 
further  analysis  of  the  seasonal  variation  of  the  NEC 
bifurcation. 

5  Bifurcation  latitude  changes  with  depth  and  season 

For  the  model,  wc  define  the  bifurcation  point  as  a 
point  of  zero  meridional  velocity,  calculated  using  90- 
day  running  mean  data  for  meridional  velocities  aver¬ 
aged  from  the  coast  to  126°  E.  Our  results  are  only 
slightly  dependent  on  this  longitude,  with  less  than 
1/4°  change  in  bifurcation  latitudes  if  the  longitude  is 
changed  by  1°.  On  the  annual  average,  the  change  is 
much  smaller.  The  bifurcation  from  our  model  is  shown 
in  Fig.  1  as  Filled  in  circles. 

As  in  the  observations,  the  model  bifurcation  lat¬ 
itude  shifts  poleward  with  depth.  Kim  et  al.  (2004) 
also  found  this  to  be  the  case  in  their  OGCM  results 
(see  their  Fig.  4).  However,  the  change  in  bifurcation 
latitude  with  depth  is  the  largest  in  the  observations  and 
the  smallest  in  the  OGCM. 


Fig.  3  Transport  of  the  NEC 
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Fig.  4  Seasonal  model  currenis  in  lhe  upper  layer  on  January  IS 
(a)  and  July  18  (b).  A  90-dav  running  average  was  applied.  The 
length  of  each  vector  is  proportional  lo  lhe  square  rool  of  the 
magnitude  of  the  currenl.  tJnii  is*  (m/s)? 


The  seasonal  cycle  of  the  bifurcation  is  mainly  con¬ 
trolled  by  a  low  dynamic  height  off  the  Philippines  coast 
near  14°  N  during  the  northern  winter  which  is  associ¬ 
ated  with  a  southward  flow  anomaly  along  the  coast. 
Figure  5  shows  the  sum  of  layer  1  and  layer  2  thick¬ 
nesses  which  is  closely  related  to  the  dynamic  height. 
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During  the  summer,  this  low  is  replaeed  by  a  positive 
dynamie  height  anomaly  and  associated  poleward  flow 
along  the  eoast,  showing  a  remarkable  agreement  with 
observations  (Qu  and  Lukas  2003).  This  annual  pres¬ 
sure  anomaly  is  caused  by  a  combination  of  remotely 
foreed  Rossby  waves  as  well  as  a  response  to  local  winds 
(Qiu  and  Lukas  1996;  Kim  et  al.  2004). 

6  Wind  stress  climatologies 

The  4.5-layer  model  was  forced  with  other  wind  stress 
climatologies  to  investigate  if  the  results  above  were 
sensitive  to  the  wind  stress  product.  As  alternative  forc¬ 
ings,  we  chose  the  Hellerman  and  Rosenstein  (1983) 
wind  stress  (HR  forcing),  a  monthly  wind  climatology 
using  QuikSCAT  winds  from  September  1999  to  March 
2002  (QS  forcing),  and  wind  stress  from  the  FSU  (1970- 
1999).  They  represent  a  range  of  wind  stress  produets 
whieh  have  well-known  advantages  and  problems.  For 
instance,  HR  forcing  is  considered  to  be  too  strong  in 
the  tropics,  but  has  been  used  in  numerous  model  in¬ 
tegrations.  QuikSCAT  scatterometer  winds  have  high 
spatial  resolution  (0.5°),  but  cover  too  few  years  to  be  a 
representative  climatology*  It  has  been  included  here 
because  of  its  high  spatial  resolution.  The  ECMWF 
ERA-15  wind  stress  used  in  the  control  is  a  relative 
coarse  (2.5°)  standard  global  reanalysis  produet,  while 
the  FSU  wind  stress  is  a  highly  quality  controlled  prod¬ 
uct  based  on  ship  observations.  Both  HR  and  FSU 
winds  have  a  spatial  resolution  of  2.0°. 

Figure  6  shows  the  Sverdrup  transport  streamfune- 
tion  calculated  for  the  four  wind  stress  products.  Note 
that  the  Hellerman-Rosenstein  wind  stress  results  in 
the  largest  transport,  followed  in  strength  by  the 
ECMWF  and  the  FSU  wind  stress  climatologies.  The 
QuikSCAT  wind  stress  is  significantly  weaker  than  the 
other  wind  stress  products  south  of  12°  N.  l  he  annual 
average  line  of  zero  Sverdrup  transport  is  at  14°  N 
for  the  QuikSCAT  wind  stress  while  it  is  about  1° 
further  poleward  for  the  other  wind  stress  products. 
The  QuikSCAT  wind  stress  used  here  is  dominated 
by  strong  La  Nina  conditions,  which  tends  to  shift  the 
bifurcation  equatorward  (Kim  ct  al.  2004). 

7  Sensitivity  to  wind  climatology 

Figure  7  shows  the  bifurcation  latitude  for  the  upper 
three  model  layers,  using  four  different  wind  stress  cli¬ 
matologies:  the  elassic  and  commonly  used  Hellerman 
and  Rosenstein  (1983)  in  Fig.  7a,  the  ECMWF  forcing 
used  as  our  foreing  for  eontrol  runs  in  Fig.  7b,  the 
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Fig.  5  Seasonal  anomaly  of 
the  sum  of  layer  1  and  layer  2 
thickness  (meters)  and 
anomaly  transports  (vectors) 
in  those  two  layers  for  DJF 
(a)  and  JJ  A  (b).  The  vectors 
has  been  rc-scaled  by  taking 
the  square  root  of  the 
magnitude  so  the  unit  is 
(nr/s)1'2 
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Fig.  6  Sverdrup  transport  for  four  different  wind  stress  products  used  in  the  simulations.  Unit  is  Sverdrup 


FSU  1979-1999  forcing  in  Fig.  7c,  and  the  QuikSCAT 
forcing  in  Fig.  7d.  The  panels  are  arranged  based  on 
the  extreme  southward  latitude  of  the  surface  layer 
bifurcation  for  a  given  climatology.  The  bifurcation 
latitude  for  each  layer  is  determined  by  taking  the  aver¬ 
age  meridional  transport  from  the  coast  to  126°  E  and 
searching  for  a  sign  change  between  10°  N  and  19°  N. 


In  layer  4,  which  represents  flow  at  about  850  m,  the 
eddy  field  dominates  the  flow,  and  a  unique  bifurcation 
is  most  often  not  identifiable  during  the  summer,  as 
also  noted  by  Qu  and  Lukas  (2003).  During  winter, 
bifurcations  in  layer  4,  north  of  the  layer  3  bifurcations, 
were  found  in  most  years,  but  variability  from  year  to 
year  does  not  provide  confidence  in  their  significance. 


Fig.  7  Bifurcation  latitude  by 
month  for  layers  1  to  3  when 
forced  by  different 
climatologies.  Layer  1:  full, 
layer  2:  dashed ,  layer  3: 
dash -dor.  From  left  to  right : 
Hcllcrman  Roscnstcin  (a), 
ECMWF  (h),  FSU  (c),  and 
QuikSC AT  (d) 


HRYR20  ECMWF  YR  20 


FSU  70  99  YR  20  QUlKscalYR20 


^  Springer 


Ocean  Dynamics  (2011)  61:1329-1344 


1337 


The  following  discussion  will  compare  the  bifurca¬ 
tion  for  the  other  three  cases  to  the  bifurcation  in  the 
control  case  forced  by  ECMWF  winds  (Fig.  7b).  For  the 
surface  layer,  the  southernmost  latitude  is  found  from 
April  through  August  and  the  northernmost  position 
in  December.  The  bifurcation  of  the  subsurface  layers 
is  closest  to  the  equator  in  June  to  August  and  at  its 
most  poleward  position  in  November  and  December. 
This  seasonal  cycle  corroborates  the  analysis  of  Qu  and 
Lukas  (2003).  Agreement  between  their  observations 
and  our  model  varies  with  depth:  It  is  at  the  observed 
latitude  in  the  upper  layer,  w  hich  is  centered  at  a  depth 
around  50  m,  and  directly  affected  by  the  wind  stress.  In 
the  second  model  layer,  which  is  mainly  geostrophieallv 
driven,  the  bifurcation  in  the  model  is  about  1°  closer 
to  the  equator  than  observed.  At  intermediate  depth, 
about  450  m,  the  observed  bifurcation  is  2°  further 
poleward  in  the  summer  and  3°  further  poleward  in 
the  winter.  The  change  with  depth  of  the  bifurcation 
latitude  is  1.5°  in  this  model  run,  but  is  a  factor  of  two 
larger  in  the  observations  (Fig.  I). 

Figure  7a  shows  that  the  amplitude  of  the  seasonal 
cycle  is  much  reduced  using  the  HR  forcing,  although 
the  phase  in  the  deepest  layer  is  similar  with  that  for  the 
ECMWF  winds.  The  equatorward  displacement  during 
the  summer  is  absent  for  the  upper  layers.  The  south¬ 
ernmost  bifurcation  takes  place  in  March  for  the  upper¬ 
most  layer  and  in  February  for  layer  2.  A  poleward  shift 
in  latitude  with  depth  is  reproduced,  but  is  less  than  in 
the  control  case.  In  contrast,  the  OS  forcing  produces 
a  model  response  with  larger  change  of  the  bifurcation 
latitude  with  depth,  about  3.5°,  and  a  bifurcation  1.5° 
closer  to  the  equator.  The  seasonal  cycle  has  a  larger 
amplitude  with  its  southernmost  position  in  June  and 
northernmost  positions  in  March  and  in  November. 
(Fig.  7d). 

The  response  to  the  FSU  forcing  (Fig.  7c)  is  in  be¬ 
tween  the  ECMWF  and  QuikSCAT  responses  in  terms 
of  the  range  of  bifurcation  change  with  depth.  How¬ 
ever,  the  amplitude  of  the  annual  cycle  of  the  bifurca¬ 
tion  latitude  is  about  the  same  for  the  four  different 
cases,  about  1.5°,  while  the  observations  by  Qu  and 
Lukas  (2003)  show'  a  range  of  2.5°. 

In  Fig.  7,  the  results  are  ordered  from  left  to  right 
using  wind  stress  climatologies  with  increasing  magni¬ 
tude  of  the  Sverdrup  transport ,  i.e..  Fig.  6.  Note  that  the 
stronger  wind  stress  curl,  i.e.,  HR  and  ECMWF  forcing, 
result  in  a  bifurcation  further  to  the  north  for  layers  1 
and  2  and  with  less  depth  dependence  than  the  weaker 
wind  stress  climatologies,  c.g.,  FSIJ  and  QuikSCAT. 

A  6-year  simulation  with  the  QuikSCAT  wind  stress 
reduced  globally  by  a  factor  of  2  (not  shown),  moved 
the  bifurcation  latitude  further  southward  by  about  1° 


for  layers  1  and  2  with  no  change  in  phase.  In  layer  3, 
the  bifurcation  occurred  less  than  0.5°  further  south, 
while  a  poleward  shift  of  about  1°  was  found  in  layer 
4.  While  the  bifurcation  of  the  depth-integrated  flow 
moved  equatorward  by  0.5°  after  reducing  the  strength 
of  the  wind  stress  by  a  factor  of  2,  the  poleward  shift  in 
bifurcation  with  depth  was  increased  as  the  strength  of 
the  forcing  was  decreased.  Wc  find  that  this  tendency 
of  increased  depth  dependency  and  equator  shift  in 
bifurcation  with  weaker  wind  stress  to  be  general  for 
our  reduced  gravity  model  results. 


8  Local  and  remote  forcing 

llie  Southeast  Asian  monsoon  has  a  dominant  wind 
component  from  the  south  during  May  through 
September,  while  the  meridional  wind  component  is 
from  the  north  in  the  remainder  of  the  year.  This 
reversal  extends  eastward  to  about  140°  E. 

The  influence  of  remote  forcing  by  annual  Rossby 
waves  is  clearly  seen  from  Fig.  8  that  displays  monthly 
maps  of  the  correlation  between  sea  surface  height 
(SSH)  with  SSHiocai,  defined  as  the  area  averaged  SSH 
over  a  small  region  bounded  by  124°  E  to  126°  E  and 
14.5°  N  to  15.5°  N  near  the  bifurcation  latitude.  A  large 
positive  correlation  approaches  from  the  northeast  of 
the  bifurcation  and  is  at  20°  N  at  lag  minus  2  months 
with  a  propagation  consistent  with  a  first  baroclinic 
Rossby  wave  signal.  Positive  lags  show  coherent  SSH 
anomalies  along  the  coast  after  the  impact  of  the 
Rossby  wave.  It  is  clear  that  remote  forcing  plays  an 
essential  role  for  the  bifurcation. 

In  order  to  determine  the  relative  importance  of 
local  and  remote  forcing,  we  define  a  local  area  bor¬ 
dered  by  longitudes  120°  E  to  140°  E  and  latitudes  5° 
N  to  20°  N.  Around  the  local  area,  we  define  a  10° 
wide  buffer  zone ,  while  the  remaining  model  domain  is 
defined  as  the  remote  area .  We  made  two  integrations: 
a  local  forcing  experiment,  where  ECMWF  monthly 
winds  were  used  in  the  local  area ,  and  annual  averaged 
wind  stress  in  the  remote  area .  In  the  buffer  zone,  a 
linear  interpolation  between  monthly  wind  and  annual 
averaged  wind  was  used.  In  the  remote  forcing  exper¬ 
iment,  annual  averaged  forcing  was  used  in  the  local 
area  and  monthly  forcing  in  the  remote  area  with  linear 
interpolation  in  the  buffer  zone. 

When  the  wind  forcing  is  allowed  only  to  vary  with 
time  within  the  local  area ,  the  winter  low  pressure 
anomaly  in  the  upper  ocean  is  weakened  compared  to 
the  control  case  (not  shown),  and  the  response  to  the 
annual  march  of  the  southeast  asian  monsoon  is  much 
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FiR.8  C  orrelation  maps  of  monthly  SSH  with  the  area  averaged 
SSHiocaj,  computed  over  the  bifurcation  region  (124°  E  126°  E, 
14.5°  N  15.5°  N)  shown  by  a  black  rectangle  in  the  panel  with  zero 
lag.  Monthly  lags  arc  shown  from  negative  7  (upper  left )  where 


SSHjocai  leads,  to  a  positive  lag  of  4  months  where  SSHj,Ka|  lags 
(lower  right).  The  solution  is  from  year  20  and  was  forced  by  the 
ECMWF  wind  stress 


faster  (Fig.  9a)  than  in  the  control  case  (Fig.  7b),  i.e., 
there  is  a  phase  shift  such  that  the  extreme  bifurcation 
latitudes  occur  about  a  month  earlier.  For  instance,  as 
the  positive  wind  stress  curl  weakens  in  the  late  winter 
just  northeast  and  east  of  the  bifurcation  region  and  the 
line  of  zero  wind  stress  curl  move  southward,  the  bifur¬ 


cation  latitude  rapidly  moves  southward  in  response.  A 
fast  return  northward  takes  place  in  November  when 
the  positive  wind  stress  curl  returns.  In  contrast,  when 
only  the  wind  east  of  140°  E  has  an  annual  cycle,  the 
bifurcation  response  is  delayed  by  1  to  2  months  in 
the  upper  ocean  (Fig.  9b).  When  monthly  forcing  is 
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Fig.  9  Bifurcation  latitude  by 
month  for  layers  1  to  3  when 
forced  by  ECMWF  wind 
stress.  Layer  1:  full,  layer  2: 
dashed,  layer  3:  dash-dot. 
From  left  to  right :  Monthly 
local  winds  and  annually 
averaged  remote  winds  (a), 
monthly  remote  winds  and 
annually  averaged  loeal  winds 
(b)  ,  as  control  run,  but  with  a 
straight  Philippine  coast  (c). 
and  as  4.5  layer  model  control 
run,  but  with  bottom 
topography  included  (d) 


Local  fofdngortyYR  20 


ECMWF  remote  forcing  YR  20 


Straight  coast  YR  20 


EC  bathy  YR  20 


applied  to  the  entire  ocean,  the  time  span  where  the 
bifurcation  latitude  is  near  its  extremes  is  significantly 
longer  (Fig.  7b).  The  Mindanao  Dome  is  southeast  of 
the  bifurcation  and  is  part  of  the  bifurcation  problem. 
Tozuka  et  al.  (2002)  demonstrated  how  local  positive 


wind  stress  curl  produced  upwelling  during  spring  was 
important  for  generating  the  Mindanao  Dome,  and  its 
damping  later  in  the  year  was  caused  by  downwelling 
generated  by  remote  negative  wind  stress  curl.  Our 
results  are  consistent  with  their  findings. 


Fig.  10  Bottom  topography 
used  in  the  five-layer  model  is 
shown  near  the  Philippine 
eoast.  I  he  depth  of  the  ocean 
varies  from  1 ,500  to  5,500  m 
in  the  model.  Contour 
interval  is  150  m 


1200  m 


by  1 50  m 


5850  m 
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9  Sensitivity  to  shape  of  coast  and  bottom  topography 

In  this  section,  the  sensitivity  to  the  shape  of  the  coast 
lines  and  to  bottom  topography  is  briefly  addressed. 
The  model  runs  with  modified  coastlines  and  bottom 
topography  were  done  with  the  same  ECMWF  wind 
forcing  as  in  the  control  run.  The  model  Philippine 
coast  line  varies  from  about  127°  E  at  6°  N  to  122.5°  E 
at  18°  N.  To  determine  if  the  shape  of  the  coastline  has 
an  influence  on  the  bifurcation,  a  model  geometry  with 
a  straight  coastline  situated  at  125°  E  was  used  in  the 
model.  The  effect  of  reshaping  the  coast  is  relatively 


small  for  the  poleward  displacement  with  depth  as  seen 
by  comparing  the  bifurcation  when  a  realistic  coastline 
is  used  (Fig.  7b)  with  the  case  when  a  straight  coastline 
is  used  (Fig.  9c).  However,  the  bifurcation  latitude 
minimum  is  delayed  until  August  and  is  0.25°  further 
to  the  south,  resulting  in  a  more  pronounced  minimum 
latitude.  During  the  fall,  the  bifurcation  latitudes  of 
the  layer  1  and  layer  2  flows  are  essentially  at  the 
same  latitude.  Overall,  the  influence  of  the  shape  of  the 
coastline  is  small. 

To  determine  the  effect  of  bottom  topography,  a 
five-layer  version  of  the  model  was  used.  "This  model 


Difference  between 
the  annual  average  flow  in 
the  five-layer  model  and  the 
4.5-layer  model  when  forced 
with  ECMWF  winds  in  layer 
1  (top)  and  in  layer  4  (bottom) 
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Annual  mean  difference  between  5  layer  and  4.5  layer  model 
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is  identical  to  the  4.5-layer  model,  except  that  a  finite 
depth  with  realistic  bottom  topography  is  used  below 
1,200  m  in  the  entire  model  domain.  The  local  bottom 
topography  is  shown  in  Fig.  10.  The  numerical  method 
used  is  the  gravity  wave  retardation  method  (Jensen 
1996)  with  a  speed  up  factor  of  8  for  the  barotropic 
gravity  waves.  This  method  has  been  shown  to  give  re¬ 
alistic  flow  over  deep  topography  (Jensen  2001, 2003). 

Introduction  of  bottom  topography  has  three  effects: 
The  finite  depth  slows  down  the  baroclinic  wave  speeds 
compared  to  the  reduced  gravity  model:  Near  155°  E, 
the  first  vertical  mode  baroclinic  planetary  wave  speed 
is  0.59  m/s  for  the  4.5-layer  model  and  0.52  m/s  for  the 
five-layer  model.  Near  the  coast,  at  130°  E,  the  phase 
speeds  have  increased  to  0.73  and  0,69  m/s,  respec¬ 
tively.  Secondly,  it  adds  a  barotropic  flow  that  in  the 
absence  of  stratification  would  follow  ///)  contours, 
where  D  is  the  ocean  depth.  Stratification  limits  the 
importance  of  local  topography,  but  can  still  affect  the 
upper  layers  through  the  joint  effect  of  baroclinieity 


and  relief  (JEBAR).  Cane  et  al.  (1998)  have  given  a 
very  useful  interpretation  that  JEBAR  has  an  effect 
on  the  barotropic  flow  when  there  is  a  near-bottom 
geostrophic  flow  up  or  down  the  topographic  slopes 
(see  their  Eq.  2),  an  explanation  that  is  particularly 
applicable  to  layer  models.  Thirdly,  interfacial  friction 
between  the  deepest  layer  and  the  upper  layers,  layer 
4  in  our  case,  is  affecting  the  strength  of  the  coupling 
between  the  deep  and  the  upper  ocean. 

For  the  upper  three  layers,  the  influence  of  topog¬ 
raphy  is  small  for  the  flow  away  from  the  boundary, 
suggesting  at  first  that  bottom  topography  effects  are 
unimportant.  However,  when  the  bifurcation  is  com 
puted,  a  shift  in  phase  is  seen  (Fig.  9d).  Figure  1 1  shows 
differences  between  the  five-layer  model  and  the  4.5- 
layer  model  simulation  for  the  annual  averaged  flow  for 
year  20  in  layer  1  and  in  layer  4.  The  differences  are 
mainly  due  to  the  additional  barotropic  component  of 
the  flow  in  the  finite  depth  model,  and  consequently, 
the  differences  between  the  two  model  solutions  in 


Fig.  12  Annual  average  flow 
in  the  abyssal  layer  during 
year  20  from  the  model 
simulation  with  bottom 
topography.  The  contours 
show  the  2,000-  and  3,000-m 
isobath  (dashed  lines)  and  the 
4,000-m  isobath  (full  line) 
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layer  2  and  layer  3  (not  shown)  are  very  close  to  those 
found  in  layer  1.  There  is  an  enhanced  westward  flow  in 
the  northern  part  of  the  NEC  and  a  decrease  in  current 
velocity  in  its  southern  part.  This  flow  anomaly  is  more 
enhanced  during  the  winter  and  leads  to  a  northward 
shift  in  the  bifurcation  latitude  during  the  fall  and 
winter  (Fig.  9d)  compared  to  the  control  ease  (Fig.  7b). 
Another  difference  between  the  finite  depth  run  and  re¬ 
duced  layer  model  run  is  in  the  phase  speed  of  westward 
propagating  Rossby  waves.  The  differences  are  fairly 
small,  and  with  the  slower  propagation  in  the  case  with 
topography,  a  lag  in  phase  of  the  bifurcation  compared 
to  the  4.5-layer  case  would  be  expected.  Rather,  the 
extrema  occur  about  a  month  earlier,  suggesting  that 
the  propagation  speed  differences  are  not  important. 
Finally,  the  topography  close  to  the  Philippine  coast 
features  a  seamount  centered  at  125°  E,  16.5°  N,  which 
has  an  anti-cyclonic  circulation  throughout  the  year. 
Figure  12  shows  this  anti-cyclonic  circulation  and  the 
abyssal  (layer  5)  annual  mean  flow  in  the  domain  of 
interest.  The  seamount  does  not  affect  the  bifurcation 
significantly,  and  as  indicated  by  the  flow  in  layer 
4  (Fig.  11,  bottom  panel),  there  is  no  indication  the 
seamount  has  any  influence  on  the  upper  layers.  The 
abyssal  flow  is  southward  along  the  Philippine  coast 
and  forms  an  anticlockwise  circulation  with  northward 
flow  along  the  western  side  of  the  Palau  Ridge  and 
westward  flow  beneath  the  NEC  The  flow  is  primarily 


along  the  isobaths  in  the  annual  mean,  but  up-  and 
downslope  flows  occur  throughout  the  year  suggesting 
some  influence  of  JEBAR.  However,  to  a  large  extent, 
the  deep  flow  is  decoupled  from  the  flow  in  the  upper 
layers. 

10  Discussion 

Earlier  modeling  work  using  wind-driven  models  (e.g., 
Qiu  and  Lukas  1996)  did  not  explore  the  change  of 
the  bifurcation  with  depth  and  found  a  seasonal  change 
that  does  not  agree  with  recent  analysis  of  observations 
(Qu  and  Lukas  2003).  As  demonstrated  in  this  paper,  a 
global  model  with  additional  layers-and  therefore  addi¬ 
tional  baroclinic  modes,  gives  results  in  good  agreement 
with  observations.  Do  the  additional  baroclinic  modes 
make  the  difference?  Experiments  with  a  1 .5-layer  ver¬ 
sion  of  the  model  reduced  the  annual  amplitude  of 
the  bifurcation  latitude,  but  did  not  change  its  phase. 
In  order  to  further  explore  this,  Miyama  et  al.  (2003) 
used  a  1.5-layer  model  with  the  same  domain  as  Qiu 
and  Lukas  (1996),  i.e..  Pacific  Ocean  only  and  identical 
forcing  (FSIJ  1961-1992;  Fig.  13).  They  found  that  the 
bifurcation  latitude  was  sensitive  to  the  model’s  domain 
si/e.  With  a  southern  boundary  at  15°  S,  a  run  was  done 
with  the  northern  boundary  at  60°  N,  and  was  compared 
with  a  run  with  the  northern  boundary  at  40°  N  the 


Fig.  13  Seasonal  variations 
of  the  bifurcation  latitude  in  a 
linear  1.5-layer  model  with 
the  northern  model  boundary 
al  60°  N  (full  line),  with 
northern  boundary  al  40°  N 
(dashed  line)  as  in  Qiu  and 
Lukas  (1996)  and  with  lhe 
northern  boundary  at  40°  N 
and  strong  damping  applied 
( dotted  line).  From  Miyama 
et  al.  (2003) 
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latitude  as  was  used  by  Qiu  and  Lukas  (1996).  The  sec¬ 
ond  run  shows  a  shift  of  the  most  equatorward  bifurca¬ 
tion  latitude  from  late  summer  to  early  spring.  The  ear¬ 
lier  response  was  due  to  a  spurious  coastal  Kelvin  wave 
that  propagated  the  effects  of  remote  winds  westward 
at  a  faster  rate  than  possible  by  Rossby  waves  (Miyama 
et  al.  2003).  They  also  made  an  experiment  with  strong 
horizontal  eddy  viscosity  along  the  boundaries  as  used 
by  Qiu  and  Lukas  (1996),  but  that  did  not  improve  the 
result  very  much  (Fig.  11).  We  therefore  suspect  the 
result  by  Qiu  and  Lukas  (1996)  was  influenced  by  their 
artificial  boundaries. 

lhe  4.5-layer  model  results  agree  quite  well  with  the 
observed  seasonal  cycle  of  the  bifurcation  latitude  and 
its  variation  with  depth.  As  in  the  analysis  by  Qu  and 
Lukas  (2003),  a  depth  of  no  motion  is  implicitly  given 
in  the  4.5-layer  model. 

Experiments  where  only  annual  mean  wind  forcing 
was  used  east  of  140°  E  provided  a  scenario  where  only 
local  seasons  had  an  influence.  This  resulted  in  shorter 
duration  of  the  time  where  the  bifurcation  latitude  was 
at  its  extremes.  The  southernmost  and  northernmost  bi¬ 
furcation  was  found  earlier  than  observed,  implicating 
that  remote  forcing  is  needed  to  get  the  correct  phase. 

An  experiment  with  straight  coastline  demonstrated 
some  impact  on  the  bifurcation,  but  indicates  that  the 
local  boundary  condition  is  of  less  importance  than 
the  forcing.  A  five-layer  model  simulation  including 
bottom  topography  produced  a  northward  shift  in  bi¬ 
furcation  during  fall  and  winter  due  to  inclusion  of 
the  barotropic  mode,  lhe  seasonal  phase  of  the  bi¬ 
furcation  was  leading  that  of  the  4.5-layer  model  by 
about  a  month,  a  result  that  is  in  less  agreement  with 
observations  than  the  4.5-layer  model  results.  However, 
inclusion  of  finite  depth  and  bottom  topography  did  not 
substantially  alter  our  results  for  the  upper  layer  flow. 

All  model  simulations  were  run  without  seasonal 
buoyancy  forcing,  but  only  a  vertical  hydromixing  that 
ensures  the  model  stratification  does  not  change  much 
from  the  initial  conditions.  This  suggests  that  the  bi¬ 
furcation  latitude  is  determined  by  the  wind  stress  and 
transferred  to  the  deeper  ocean  via  baroelinie  waves. 
Different  wind  climatologies  produce  qualitatively  sim¬ 
ilar  results,  but  also  show  one  important  difference: 
Stronger  wind  stress  results  in  a  bifurcation  further 
poleward. 

Finally,  it  should  be  noted  that  inter-annual  variabil¬ 
ity  of  the  bifurcation  latitude  associated  with  ENSO  has 
been  found  in  modeling  studies  and  has  recently  been 
observed.  Modeling  studies  by  Kim  et  al.  (2004)  found 
a  northward  shift  in  the  position  of  the  bifurcation 
during  El  Nino  and  a  southward  shift  during  Lit  Nina. 
Observations  by  Kashino  et  al.  (2(X)9)  did  not  find  a 


significant  difference  in  the  bifurcation  latitude  during 
the  2006-2007  El  Nino  and  the  2007-2008  La  Nifta, 
but  did  find  increased  North  Equatorial  Current  and 
Mindanao  Current  transports. 

Overall,  the  seasonal  cycle  of  the  bifurcation  latitude 
is  fairly  robust,  and  the  annual  average  bifurcation 
latitude  differs  less  than  0.5°  between  the  runs  forced 
with  different  climatological  wind  stress  products.  An 
exception  is  the  run  forced  with  QuikSCAT  winds. 
However,  this  wind  forcing  was  chosen  as  a  short  non¬ 
representative  climatology,  suggesting  that  the  bifurca¬ 
tion  is  subject  to  inter-annual  variability. 
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